On the theory of superconductivity in the extended Hubbard model: 

Spin-fluctuation pairing 
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A microscopic theory of superconductivity in the extended Hubbard model which takes into 
account the intersite Coulomb repulsion and electron-phonon interaction is developed in the limit 
of strong correlations. The Dyson equation for normal and pair Green functions expressed in terms 
of the Hubbard operators is derived. The self-energy is obtained in the noncrossing approximation. 
In the normal state, antiferromagnetic short-range correlations result in the electronic spectrum 
with a narrow bandwidth. We calculate superconducting T c by taking into account the pairing 
mediated by charge and spin fluctuations and phonons. We found the d-wave pairing with high- 
T c mediated by spin fluctuations induced by the strong kinematic interaction for the Hubbard 
operators. Contributions to the d-wave pairing coming from the intersite Coulomb repulsion and 
phonons turned out to be small. 

PACS numbers: 74.20.Mn, 71.27. -fa, 71.10.Fd, 74.72.-h 



I. INTRODUCTION 

Despite intensive studies of high-temperature super- 
conductivity (HTSC) in cuprates for many years af- 
ter the discovery of Bednorz and Miiller [l|, a com- 
monly accepted mechanism of HTSC is still lacking (see, 
e.g. @, Q). A good candidate from various proposed 
mechanisms is based on a model of strongly correlated 
electrons [1J. In the model, superconductivity occurs at 
finite doping in the resonating valence bond state (RVB) 
due to the antiferromagnetic (AF) superexchange in the 
t-J model. A possibility of HTSC mediated by AF spin 
fluctuations as a "glue" for superconducting pairing was 
also considered [f| , mostly within phenomenological spin- 
fcrmion models (see, e.g., [6|-|9(, and references therein). 

Recent studies of spin-excitations by magnetic inelas- 
tic neutron scattering (INS) and the electronic spectrum 
py angle-resolved photoemission spectroscopy (ARPES) 
have revealed an important role of AF spin excitations 
in the "kink" phenomenon and the d-wave pairing in 
cuprates (see, e.g. , [To| and references therein). In par- 
ticular, in Ref. using INS and ARPES studies on the 
same YBa2Cu3C>6.6 (YBCOg.6) crystal, an estimation for 
superconducting T c ~ 150 K was found. The main argu- 
ment against the spin-fluctuation pairing, the weak inten- 
sity of spin fluctuations at the optimal doping seen in INS 
experiments [l2| , was dismissed in the recent resonant in- 
elastic x-ray scattering experiments [131 ] . In a large family 
of cuprate superconductors, paramagnon AF excitations 
with the dispersion and spectral weight similar to those 
of magnons in undoped cuprates were observed. Using 
the magnetic spectrum found in the YBCO7 crystal, su- 
perconductivity with T c = 100 — 200 K was predicted. 
Thus, spin fluctuations have sufficient strength to me- 
diate HTSC in cuprates and to explain various physical 
properties of cuprate materials as, e.g., the optical con- 
ductivity [lij • Therefore, it can be suggested that the al- 



ternative mechanism based on the conventional electron- 
phonon interaction (EPI) (see, e.g., [HI, EH) plays a sec- 
ondary role in the cuprate superconductors. 

Recently, in Ref. [l7j using the renormalization group 
(RG) method an asymptotically exact solution for the 
d-wave pairing was found in the conventional Hubbard 
model [18] in the weak correlation limit, U -C t. How- 
ever, as was pointed out later in Ref. [19}], a contribution 
from the repulsive well-screened weak Coulomb interac- 
tion (CI) in the first order strongly suppresses the pairing 
induced by the contributions of higher orders, and a pos- 
sibility for superconductivity "from repulsion" was ques- 
tioned. At the same time, in Ref. [2(| it was shown that 
the p-wave superconductivity exists in the electronic gas 
at low density with a strong repulsion U and a relatively 
strong Coulomb intersite interaction Vjj (see, also [21| 
and references therein) . Later on, in Ref. [22J RG studies 
of the extended Hubbard model with the intersite inter- 
action have shown that superconducting pairing of vari- 
ous symmetries, extended s-, p-, and d-wave types, can 
occur depending on the electron concentration and the 
intersite interaction Vij. However, in these investigations 
the Fermi-liquid model in the weak correlation limit was 
used. To study superconductivity in cuprates, the Mott- 
Hubbard (more accurately, charge-transfer) doped insu- 
lators, a theory of strongly correlated electronic systems 
should be used (for reviews see [H, H3])- 

In the present paper we consider superconductiv- 
ity in the extended Hubbard model with a weak in- 
tersite Coulomb repulsion Vij but in comparison with 
Refs. [13, EH HHh we study the limit of strong corre- 
lations, U S> t. To compare various contributions to 
the superconducting d-wave pairing, we consider also a 
model of the EPI with strong forward scattering proposed 
in Ref. [25]. The Dyson equation for the thermodynamic 
Green functions (GFs) expressed in terms of the Hubbard 
operators (HOs) is derived using the Mori-type projec- 
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tion technique 26]. The self-energy is calculated in the 
noncrossing approximation (NCA) as in the microscopic 
theory of the electronic spectrum in the normal state in 
our previous publication [2?} ■ We show that the kine- 
matic interaction for the HOs generates the AF superex- 
change pairing similar to the t— J model. A contribution 
from the intersite Coulomb repulsion in the first order 
suppresses the pairing as found in Refs. But 
the kinematic interaction induces also a strong electron 
interaction with spin-fluctuations which results in the d- 
wave superconductivity with high-T c . Contribution from 
the EPI to the d-wave pairing turned out to be small. 

In the next section we introduce the model, derive 
the Dyson equation, and calculate the self-energy in the 
NCA. A self-consistent system of equations is formulated 
in Sec. Mil Results of computations of the electronic spec- 
trum in the normal state and of superconducting T c and 
the d-wave gap function are presented in Sec. IIVI Con- 
cluding remarks are given in Sec. [V] In the Appendix 
details of the calculations are given. 



II. GENERAL FORMULATION 
A. Extended Hubbard model 

We consider an extended Hubbard model on a square 
lattice which we write in terms of the HOs [2811: 
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To apply the model for consideration of the cuprate su- 
perconductors, we introduce the HOs for holes taking 
into account four possible states on a lattice site i: an 
empty state (a,/3 = 0), a singly occupied hole state 
(a, f3 — a) with the spin a/2 = ±(1/2) , a = —a, and a 
two-hole state (a,/9 = 2). Then the HO xf = \ia){i0\ 
describes the transition from the state \i,/3) to the state 
\i,a). Energy parameters in the model (JlJ are taken 
close to the values found within the cell-perturbation 
method [29| for the p-d model for the Cu02 plane [30| . 
In particular, the single-particle energy e± = £d — H is 
the energy of the <i-type one-hole state measured from 
the chemical potential fi and the two-particle energy 
82 = 2ei + U is the energy of the two-hole p-d Zhang- Rice 
singlet state [3l[ . The effective Hubbard U in cuprates is 
the charge-transfer energy U = A p d = e p ~ Accord- 
ing to the cell perturbation method, in general case the 
values of the hopping parameters t^j in (fT|) depends on 
the the subband indices n, £ = 1, 2. 

In the last term in (fTJ) in addition to the inertsite CI 
Vij for holes in the plane we take into account also the 
EPI gij for holes with phonons: 



Hc,e P = \ E VijNiN) + ]T 



9ij 



N, 



j ■■ 



(2) 



where Uj describes an atomic displacement on the lat- 
tice site j for a particular phonon mode. More generally, 
the EPI can be written as a sum ^ y gfjUj over the nor- 
mal modes v. The hole number operator and the spin 
operators in terms of HOs are defined as 

Ni = 2^^r + 2Xf, (3) 

cr 

Sf = X[-, St - (a/2) [Xr - X?% (4) 

The completeness relation for the HOs, X®° + X° a + 
Xf = 1, rigorously preserves the constraint of no double 
occupancy of the quantum state a on any lattice site i. 
From the multiplication rule X^X? 5 = 5p 1 X" 5 follows 
the commutation relations: 



y Xf ± 8 5a xf 



(5) 



The upper sign refers to the Fermi-type operators like 
X° cr while the lower sign refers to the Bose-type operators 
like N{ (]3]) or the spin operators (jU). 

The chemical potential /i depends on the average hole 
occupation number 



1 



'.Ni. 



(G) 



where (...) denotes the statistical average with the Hamil- 
tonian (fTJ). 

We emphasize here that the Hubbard model flU does 
not involve a dynamical coupling of electrons (holes) to 
fluctuations of spins or charges. Its role is played by the 
kinematic interaction caused by the complicated com- 
mutation relations ([5]), as was already noted by Hub- 
bard [28|. For example, the equation of motion for the 
HO Af 2 in the Heisenberg representation (H = 1) reads, 



% -X^ 
di 1 



= [Xf,H] = (e 1 + U)Xf 

+ Y.i^BH^-a^Bl^xf 

l,a' 



i 



Here BV~ , are the Bose-like operators, 



D 



21 



{X? + xr)5o'a+Xf5 a , s 
{N l /2 + aS z l )8 a , CJ + S° 5 a , s , 
{Ni/2 + aSf)S a/(T - S?6 a , 9 , 



(J) 



(8) 



(9) 



where we used the definition of the number operator ^) 
and the spin operators (0} . The last term in (|7J| is caused 
by the dynamic intersite CI and the EPI. 



B. Dyson equation 

To consider the superconducting pairing in the model 
(TTJ), we introduce the two-time thermodynamic GF [32| 
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G°(k,w) = (W -E CT (k)) V 



expressed in terms of the four-component Nambu opera- and the corresponding zero-order GF 
tors, X tlJ and tj a = Xf° Xf X^) : 

Ga^t-H) = -i6(t-t'){{x itT (t),xl(t')}) 

= ({x icr (t)\xl(t'))}, (10) 

where {A, B} = AB+BA, A(t) = exp(iHt)Aexp(-iHt), 
and 6{x) — 1 for x > and 6(x) = for x < 0. The 
Fourier representation in (k, w)-space is defined by the 
relations: 

1 

2^ 
1 



(19) 



where fo is the 4x4 unit matrix. 

To calculate the multiparticle GF {{zf*\t) | Xj^t'))) 
in (|16|) we differentiate it with respect to the second time 
t' and apply the same projection procedure as in (|T7|) . 
This results in the equation for the GF (|16l) in the form, 



J ij<7 (^ 

G 



i') = 



dte-^*- t ')Gy cr (a;) ) (11) 

to 

exp[ik(i-j)]G CT (k,o;). (12) 



Go - (k, cj) = G°(k,w) + G°(k,w) T CT (k, w ) G°(k, 



where the scattering matrix 



The GF (1111) is convenient to write in the matrix form 



G ija (uj) 



Gija (w) 



Gjia(-w) 



(13) 



where the normal Gij a and anomalous (pair) Fij a GFs 



are the 2x2 matrices for two Hubbard subbands: 
Xf 1 



xr 
xr 



via tactOh 



(14) 



(15) 



To calculate the GF (|10[) we use the equation of motion 
method. Differentiating the GF with respect to time t, 
the Fourier representation of it leads to the equation 



where Q — 



uG ij(7 (oj) = <%Q + (([X i(T ,H] | X}J U . (16) 

({X ia ,xl})=f xQ 
and fo is the 2x2 unit ma- 



Here the correlation function Q 
Q 2 
Qi 

trix. The spectral weights of the Hubbard subbands in 
the paramagnetic state Q2 = (X 22 + X[ a ) = n/2 and 
Qi = (Xf° + Xf s ) =1-Q 2 depend on the occupation 
number of holes © . In the Q matrix we neglect anoma- 
lous averages of the type (X® 2 ) which are irrelevant for 
the d-wave pairing [331 ]. 

To introduce the zero-order quasiparticle (QP) excita- 
tion energy we use the Mori- type projection method [26| . 
In this approach, the many-particle operator Zi a = 
[Xi a , H] in (|16p is written as a sum of a linear part and 

an irreducible part Z\ a ' orthogonal to X\ a : 



Z ia = [X i(J ,H] 



E 

1 



^ilaXi a + Z 



The orthogonality condition ({Z^\ Xj a }) = deter- 
mines the excitation energy in the mean-field approxi- 
mation (MFA) 



(ir) 



(17) 



= ({[X ia ,H},xl DQ" 1 



(1/JV) ^exp[ik(i-j)]E CT (k), (18) 

k 



Q- 



>(ir) 1 ^(ir)t 



k<r 



Q" 



(20) 



(21) 



Now we can introduce the self-energy operator E (T (q, uS) 
as the proper part (pp) of the scattering matrix (|2"Tj) 
which has no parts connected by the zero-order GF (fl9|) 
according to the equation: T = T. + T. G° T. The defini- 
tion of the proper part of the scattering matrix (|21[) is 
equivalent to an introduction of a projected Liouvillian 
superoperator for the memory function in the conven- 
tional Mori technique (26|. 

Using the self-energy operator instead of the scattering 
matrix in Eq. (|20[) we obtain the Dyson equation for the 
GF P): 



G ff (k,w) 



[UJTq 



E cr (k)-Qi: ff (k,w)]- 1 Q J 



where the self-energy operator is given by 



QE <r (k, W ) = «^|i ] 



K ir )t\\(pp) 



Q- 



(22) 



(23) 



Dyson equation (f22|) with the zero-order QP excitation 
energy (TTg)) and the self-energy (|2"3"|) gives an exact repre- 
sentation for the GF (Ti7j|) . The self-energy takes into ac- 
count processes of inelastic scattering of electrons (holes) 
on spin and charge fluctuations due to the kinematic in- 
teraction and the dynamic intersite CI and the EPI (see 
Eq. (O). To obtain a closed system of equations, the 
multiparticle GF in the self-energy operator (f2"3"|) should 
be evaluated as discussed below. 



III. SELF-CONSISTENT SYSTEM OF 
EQUATIONS 

A. Mean-field approximation 

The superconducting pairing in the Hubbard model 
already occurs in the MFA and is caused by the kinetic 
superexchange interaction as in the t-J model [|| . There- 
fore, it is reasonable to consider at first the MFA de- 
scribed by the zero-order GF ([T9"l) . Using commutation 
relations ([5]) for the HOs we calculate the energy matrix 
(HI: 



-%ja 



A 



(24) 
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The matrix e(k) = J^j exp[ik(i — j)] iij after diagonal- 
ization determines the QP spectrum in the two Hubbard 
subbands in the normal state (for details see [13] ) : 

ei )2 (k) = (l/2)[^ 2 (k)+ Wl (k)] T (l/2)A(k), (25) 
w t (k) = 4ta t7 (k) + 4A*Y(k) + 4 ( t i / Y / (k) 

+ 6^(^ + 17^,2- A*, (t=l,2) (26) 
A(k) = {[ W2 (k)- Wl (k)] 2 +4^(k) 2 } 1 / 2 , 
W(k) = 4ta 12 7(k)+4t'/3 127 '(k)+4i"/3 127 "(k). 

Here the hopping parameters in (JTJ) are assumed to be 



equal, tff = t 
expression: 



fii 



f!2 



Uj, where Uj is defined by the 



Uj = (1/N) ^exp[ik(i-j)]i(k), (27) 

k 

i(k) = 4t 7 (k) +4t' 7 '(k) +4i" 7 "(k). (28) 

The hopping parameters are equal to t for the near- 
est neighbor (n.n.) sites a\ = (±a x ,±a y ), t' 
for the next nearest neighbor (n.n.n.) sites ad = 
±(a x ± a y ), and t" - for the n.n.n. sites a 2 = 
±2a x ,±2a y . The corresponding k-dependent functions 
are: 7 (k) = (l/2)(cosfc;E + cosfc a ), 7 '(k) = cos k x cos k y , 
and 7 "(k) = (l/2)(cos2fc 2: +cos2fc a ) (the lattice con- 
stants a x ~ a y are put to unity). The contribution from 
the CI Vij in ([2"6"]l is given by 



uj 



(c) 
1(2) 



(k) 



q)^i( 2) (q), 



(29) 



where A^q) = (X°*X™ )/Q u N 2 (q) = (X- 2 X 2 -)/Q 2 
and F(k — q) is the Fourier transform of Vij. 

The kinematic interaction for the HOs results in 
renormalization of the spectrum (|25[) determined by 
the parameters: a L = Q L [l + Ci/Q 2 ], /3 L = Q L [1 + 
C 2 /Q 2 } , a 12 = VQI&[l-Ci/QiQ 2 ], (3 12 = 
C 2 /Q\Q 2 ] . Here we take into account the renormaliza- 
tion caused by the spin correlation functions for the n.n. 
and the n.n.n. sites, respectively: 



Ci = (SiS 



i+ai i 



C 2 = (SiS 



i J i+a d / 



(30) 



The short-range AF correlations considerably suppress 
the n.n. hopping parameters since C± < and at low 
doping \C\\ = 0.1 — 0.2 that results in ct t < 1. At the 
same time, the n.n.n. hopping parameters are increased 
since C 2 > 0. 

Now we evaluate the anomalous components A,-j CT of 
the matrix (|24[) which determine the superconducting gap 
in the MFA. Considering the singlet (i-wave pairing, we 
calculate the intersite pair correlation functions. The di- 
agonal matrix components are given by the equations: 



22 / 

11 C 



= -at 2 }{Xf 2 Nj)-V i j{XfX? 



at} 2 (NjX^)-Vij(X, 



os X° a ) 



(31) 
(32) 



Here we used the original notation for the interband 
hopping parameters t\ 2 to emphasize that the kinematic 



pairing (X\ Nj) is mediated by the interband hopping. 
In terms of the Fermi operators = X® a + <rX° 2 , 
the pair correlation function in pip can be written as 
(XfN 3 ) = (X° l Xf 2 Nj) = ( ail a lt N 3 ). This representa- 
tion shows that the kinematic pairing occurs on a single 
lattice site but in two Hubbard subbands (34] ]. 

The correlation function (X® 2 Nj) can be calculated di- 
rectly from the GF Lijit-t') = ((AT° 2 (i) | Nj(t'))) with- 
out any decoupling approximation as shown in Ref. [34j . 
In particular, under hole doping, n = 1 + 8 > 1, the pair 
correlation function in the two-site approximation reads 
(for details see Appendix A): 



(X° 2 N 3 ) - f . 



4t 12 

~TT a ' * X 



2 V(T2\ 



(33) 



As a result, the equation for the superconducting gap in 
pTjl can be written as 



A 22 a = (J l3 -V t3 )(X 



a2 Xf)/Q 2 



(34) 



where = 4 (t\ 2 ) 2 /U is the AF superexchange interac- 
tion. A similar equation holds for the gap in the one-hole 
subband: Ajj a = (J - - Vij) (Xf Xf) /Qx . We thus 
conclude that the pairing in the Hubbard model in the 
MFA is similar to the superconductivity in the t-J model 
mediated by the AF superexchange interaction J,j. 



B. Self-energy operator 

Self-energy operator (|23[) can be written in the same 
matrix form as the GF (IT3l): 



QZi ja (u) = I . f . Q , (35) 

where the matrices M and $ denote the respective nor- 
mal and anomalous (pair) components of the self-energy 
operator. The system of equations for the (4 x 4) ma- 
trix GF (fT3|) and the self-energy l[35|) can be reduced to a 
system of equations for the normal G,j(k, uj) and the pair 
F CT (k, uj) (2 x 2) matrix components. Using representa- 
tions for the energy matrix (|24|) and the self-energy (l35l , 
we derive for these components the following system of 
matrix equations: 

G(k,uj) = (G^k,^)- 1 

+ 6 r 0cw)G A r(k,-w)ft(k,a;)) * Q, (36) 
F a (k,uj) = -G N (k,-uj)<p a {k,uj)G{k,uj), (37) 
where we introduced the normal state GF 

G N (k,uj) = (wTb-e(k)-M(k,w)/0)~\ (38) 
and the superconducting gap function 

^(k,w) = A CT (k) + $ CT (k,w)/Q. (39) 
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To calculate the self-energy matrix (|35|) we use the mode- 
coupling approximation which is equivalent to the NCA 
in the diagram technique for GFs. In this approximation, 
a propagation of Fermi-like excitations described by the 
operator Xf 2 , and Bose-like excitations described by 
the operator Bi aa i for I ^ i is assumed to be indepen- 
dent. Therefore, the time-dependent multiparticle corre- 
lation functions in the self-energy operators (|35l) can be 
written as a product of fermionic and bosonic correlation 
functions, 



(Xf/'B^ „\Bi, 



■(t)xf 2 (t)) 



= 8„ w , (X 2 /xf 2 (t)) {B) aa , \B iaa , (t)), 

' /2d " |D >(t)xf 2 (t)) 



(Xfi Bjaa"\Bi, 



6 a/ ^(Xpxf 2 (t)) {B ]rTS ,B iaa ,{t)) 



(40) 
(41) 



The time-dependent correlation functions are calculated 
self-consistently using the corresponding GFs (for details 
see Appendix B). 

In particular, the normal and anomalous diagonal com- 
ponents of the self-energy for the two-hole subband are 
determined by the expressions 

+oo 

M 22 (k, W ) = 1]T /d*tf<+>(w,z|q,k-q) 

1 -oo 

x |-ilm[G 22 (q,z) + G 11 (q,z)]|, (42) 



+ 00 



d> 22 (k,w) = 1]T /dz *|q,k-q) 

1 -oo 

x {-~Im[F™(< l ,z)-F?(< l ,z)]y (43) 

where G a "(q, z) and F" a (q, z) are given by the diagonal 
components of the matrices (|3"6")l . (|3"T)) . Similar expres- 
sions hold for the self-energy components M 11 (k, of) and 
$3. (k, w) for electron doping when the Fermi energy lo- 
cated in the one-hole subband (see Ref. [35j]). Note, that 
in the paramagnetic normal state the GF (f31)|) and the 
self-energy (|42l do not depend on the spin a. 

The kernel of the integral equations (|4"2"j) . (|4"3")l has a 
form, similar to the strong-coupling Migdal-Eliashberg 
theory [MII3: 



K (±) (w,z|q,k- q) = - J dn 



1 + N(Q) -n(z) 
uj — z — Q 



: {|t(q)| 2 Im Xs /(k - q, A) ± |<?(k - q}\ 2 Im Xph (k - q, ft) 



±[|V(k- 



|i(q)| 2 /4]lm Xc/ (k-q,ft)}, (44) 



where »(w) = [e w / T + l]" 1 and N(uj) = [e u ' T - l]" 1 . 
The spectral densities of bosonic excitations are deter- 



mined by the dynamic susceptibility for spin (s/), num- 
ber (charge) (c/), and lattice (phonon) (p/i) fluctuations 



Xsf(q.,u) 
Xc/(q,w) 

x P ?i(q,w) 



-((S q |S_ q )) w , 

-«<W q |«5iV_ q )) 



-((Uq|U-q» W , 



(45) 

(46) 
(47) 



iV q — (iV q ), and lattice 



which are defined in terms of the commutator GFs [32 
for the spin S q , number SNq 
displacement (phonon) u q operators. 

In the NCA, vertex corrections are neglected as in the 
Migdal-Eliashberg theory. For the EPI g(k — q ) th e ver- 
tex corrections are small, as shown by Migdal [3Q|. The 
interaction i(q) with spin-fluctuations (|45)) induced by 
the intraband hopping is not small and vertex corrections 
may be important. However, as was shown in Ref. [38| 
a certain set of diagrams, in particular the first crossing 
diagram, vanishes due to kinematic restrictions for spin 
scattering processes. Moreover, in Ref. [39| it was found 
that the renormalization of the vertex for a short AF cor- 
relation length is weak. Therefore, the NCA for the self- 
energy calculated self-consistently can be considered as 
a reasonable approximation. This approximation makes 
it possible to consider the strong coupling regime which 
is essential in study of renormalization of quasiparticle 
spectra and the superconducting pairing. 



C. Two-subband model 

In this section we derive a self-consistent system of 
equations for the GFs (j36j) -(j38 l) and the self-energy com- 
ponents (|42[) . (1431) for the two Hubbard subbands adopt- 
ing several approximations to make the system of equa- 
tions numerically tractable. 

At first we consider equations for the normal state. 
The diagonal components of the GF ([38]) can be written 
as [H: 



G# (2a) (k>w) = [l-6(k)]G 1(2) (k,w) 
+ 6(k)G 2(1) (k,w), 
1 

w - £ i(2)(k) - E(k,u;) 



Gi( 2 )(k,w) = 



(48) 
(49) 



where the hybridization parameter &(k) = [^(k) — 
W2(k)]/[e2(k) — £i(k)] . The self-energy can be approxi- 
mated by the same function for the both subbands, 

+00 

E(k,u;) = /dz^+)( W ,z|q,k-q) 

1 -00 

x [-(l/7r)]Im[Gi(q,z) + G 2 (q,z)]. (50) 

The chemical potential \x is calculated from the equation 
for the average hole occupation number ([5]): 

n = 1 + 6 = 2{XD + 2(A 22 ) = A £ iV,(q), (51) 



G 



where the hole occupation number is given by 

%)(k) = %i)(k)+JV (M) (k), (52) 
N {hl) (k) = [g 1 + (n-l)6(k)]7V 1 (k), 
N m (k) = [Q 2 -(n-l)6(k)]2V 2 (k), 

N 1{2) (k) = -- y^_ 757 _l m G 1(2) (k, W ).(53) 
Density of states (DOS) is determined by 

^) = ^E A w( k >")> ( 54 ) 

k 

where the spectral function for holes reads 

A (h) (k,u) = [Q 1 +P(k)]A 1 (k,u) 

+ [Q 2 -P(k)]A 2 (k,u>), (55) 
j4i(2)(k,w) = -(l/7r)ImGi( 2 )(k,cj). 

Here the hybridization parameter P(k) = (n — l)&(k) — 
2y/Qi Q 2 [W (k) / A(k)] takes into account contributions 
from both the diagonal and off-diagonal components of 
the GF (pp. 

Now we formulate equations for the superconducting 
gap (|3"5)) . We consider the case of hole doping when the 
Fermi energy is located in the two-hole subband, n = 
1 + 5 > 1 . By taking into account the gap equation (|34|) 
in the MFA and the self-energy (|33]l . Eq. ((39j) for the 
two-hole subband gap </s(k, cj) — cr<yj 2i(7 (k, u;) reads, 

+00 

= ^E/ ^[-^Imif 2 (q,z)] 



{[J(k-q)-y(k-q)]itanh^ 

'(-) 



Jf(-)(w,z|cbk-q3}. 



(56) 



Here the contribution from the one-hole subband 
F* 1 (q, z) in (|43|) was neglected since this filled band much 
below the Fermi level gives a vanishingly small contribu- 
tion to the pairing. To determine the superconducting T c 
it is sufficient to solve a linear equation for the gap (f5B|) 
using the linear approximation for the pair GF (|37[) . 

Ff (k,w) = -G^fk, -w)G^(k,w)(7 V (k,w)Q a 
« -[l-6(q)] 2 G 2 (k,- w )G 2 (k, w ) f 7^(k, w )Q 2 . (57) 

As in (1551) . here we neglect the GF Gi(k, w) in (|4"5)l since 
the contribution to the pairing from the one-hole subband 
much below the Fermi energy is small. 

For numerical calculations, it is convenient to intro- 
duce the imaginary frequency representation, iu) n — 
inT(2n + 1), n = 0, ±1,±2, .... In this representation 
the self-energy (|5U|) reads 



S(k, uj n ) 



T 

'N 



^^A (+) (q,k-q | u n - u m ) 



x [Gi(q,w m ) + G 2 (q,w m )] 

= iw n [l-Z(k,w„)] + A:(k,w n ). (58) 



Here we introduced the renormalization parameters 

w [1 - Z(k, u)] = (1/2) [E(k, w) - S(k, -w)], (59) 
X(k, w) = (1/2) [£(k, w) + E(k, -w)]. (60) 

The normal GF (|4l?)) for the two subbands takes the form: 

{Gl(2)(k, Un)}" 1 = iu n ~ £l( 2 )(k) - E(k, w„) 

= iu; n Z(k, u n ) - [e m (k) + X(k, w„)] . (61) 

The hole occupation number (j5"3")) in terms of the GF (joTjl 

reads: 



iv 1(2) (k) = i + r^G 1(2) (k, Wm ). 



(62) 



The gap equation (|56[) in the linear approximation for 
the pair GF ([57} can be written as 

v {k,u n ) = |]T£{j(k-q)-y(k- q ) 

q m 

+ A ( ~ } (q,k- q | w n -w m )} (63) 

. [1 - fr(q)] 2 i(q,^m) 

[cj m Z(q,o; m )] 2 + [e 2 (q) + Xq,uj m )] 2 ' 

The interaction functions in (|58[) and (1631) in the imagi- 
nary frequency representation are given by 

A (±) (<1 k - q|i/ n ) = -|t(q)| 2 Xs/(k - q, v n ) 
T { [\V(k - q)| 2 + |t(q)| 2 /4] Xcf(k - q, i/ n ) 
+ |. 9 (k-q)| 2 X ^(k-q,^)}. (64) 

Thus, we have derived the self-consistent system of equa- 
tions for the normal GF (|6ip . the self-energy (|58|) . and 
the gap function (|B"3")) . 

IV. RESULTS AND DISCUSSION 

A. Model parameters 

To perform numerical calculations we should specify 
model parameters in the derived system of equations. For 
the intersite CI V(q) we consider two models. In the first 
model the CI is determined by the repulsion of two holes 
on the n.n. lattice sites, 



Vi(q) = 2Vi (cosg x + cosq y ). 



(65) 



According to the cell-perturbation method [231, f° r the 
conventional values of electronic parameters in the p-d 
model for the Cu0 2 plane, the CI of two n.n. holes is 
estimated by the value V\ = 0.1—0.2 cV. The CI for 11. n.n. 
holes, 4V 2 cos q x cos q y , is much smaller, V 2 /Vi ~ 0.04 
and can be safely neglected. 

As the second model we consider the 2D screened CI 
suggested in Ref. (l9j : 



Vc{<C 



1 



1 



ae aq 



an 



aq 



(66) 
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where a is the lattice parameter (below we put a = 1) 
and the dielectric constant eo takes into account lattice 
polarization induced by ligand fields. We assume that 
the screening parameter k depends on the doping and 
can be described by the interpolation formula: an = 4<5, 
so that an = 0.2 in the underdoped case (8 = 0.05) 
and aK = 1 for the overdoped case (<5 = 0.25). The 
energy u c we estimate from calculation of the CI (|6T>|) at 
K = for two n.n. holes at the distance a x assuming 
it to be equal to V\ in the model (|55|) : V c i(k = 0) = 
(u c /N) J2q( cos Qx / <l) — Vl ■ From this equation we get 
an estimation u c — Vi/0.175 ~ 1 eV or u c — 2.5 i for 
t = 0.4 cV. Here for convenience, we take V\ — 0.175 eV. 

In the present study we do not perform self-consistent 
computation of spin and charge excitation spectra but 
adopt certain models for the spin (|4"5"j) . charge PS) , and 
phonon (|4T|) susceptibility in Eq. (l4"4l or Eq. (f6"4"|) . Since 
we consider the electronic spectrum only in the normal 
state and calculate superconducting transition temper- 
ature T c from the linearized gap equation (|56p or (|63[) 
the feedback effects caused by opening a superconduct- 
ing gap are not essential which justifies usage of model 
functions for the susceptibility. 

Due to a large energy scale of charge fluctuations, of 
the order of several t, in comparison with the spin ex- 
citation energy of the order of J, the charge fluctuation 
contributions in Eq. (|44"]) can be considered in the static 
limit 

Xc/ (k) = Xc/,i(k) + Xc/, 2 (k), (67) 
1_ N hm (q + k) - AT ftl(2) (q) 



Xc/,l(2)( k ) 



N 



£i (2 )(q + k) - £i (2 )(q) 



where the hole occupation numbers A^ ?ll ( 2 )(q) are defined 
in Eq. (f5"2"]) . It is assumed that the system is far away 
from a charge instability or a stripe formation when the 
energy dependence of the charge susceptibility may be 
essential (see, e.g., Refs. (ioMi^ V 

For the dynamical spin susceptibility Xs/(q, <+>) 



-«s„ 

cal studies 



we take a model suggested in numeri- 



Imx s /(q, u + i0 + ) = Xs/(q) x"f(u) 

*Q tanh " 1 

l + C 2 [l + 7(q)] 2T 1 + (oj/^y 



(68) 



The q-dependence in Xs/(q) is determined by the AF 
correlation length £ (in units of a). The frequency de- 
pendence is determined by a broad spin-fluctuation spec- 
trum x'sfi 1 ^) with a cut-off energy of the order of the ex- 
change energy uj s ~ J. This type of the spin-excitation 
spectrum was found also in the microscopic theory for 
the t-J model in Ref. (44] . The strength of the spin- 
fluctuation interaction given by the static susceptibility 
Xq = Xs/(Q) at the AF wave vector Q = (n,Tr), 



XQ 



3(1-*) 
2uj s 



I N ^ 1 



1 



£ 2 [i + 7(q)] 



(69) 



TABLE I: Spin correlation functions Ci, C2, spin susceptibil- 
ity \Q, and projected spin susceptibility Xsf for several values 
of AF correlation length £ related to hole concentration 5. 





S 


Ci 


c 2 


XQ ■* 


-Xsf -t 


3.4 


0.05 


- 0.26 


0.16 


29.5 


1.32 


2.4 


0.10 


- 0.20 


0.11 


12.6 


1.05 


1.5 


0.25 


- 0.12 


0.05 


6.8 


0.61 



is fixed by the normalization condition: 

+00 



Spin correlation functions C\, C2 (|5D|) in the single- 
particle excitation spectrum (1251) can be calculated using 
the same model (l68l): 



C 



i = ^E^7(q), C 2 = l^C q7 '(q). 



(70) 



Here the spin correlation function C q = (S q S_ q ) = 

C(0/{1 + ? 2 [1 + 7(q)]} where C(0 = XQ (w./2). The 
results of computation of the correlation functions at sev- 
eral values of the AF correlation length £ related to the 
hole concentrations 8 are given in TableUwhere the static 
susceptibility xq 5 the projected spin susceptibility Xsf 
(see Eq. ([82l ) are also given. 

To estimate the contribution from phonons in Eq. (|44[) 
we consider a model susceptibility for optic phonons and 
the EPI matrix element in the form similar to Ref. 14511: 



V ep (q,u) = \g(q)\ 2 x P h(q,u) = g ep ■ 



l4 



S(q), (71) 



where g ep is the "bare" matrix element for the short- 
range EPI, while the momentum dependence of the EPI 
is determined by the vertex correction S(q) . It takes 
into account a strong suppression of charge fluctuations 
at small distances (large scattering momenta q) induced 
by electron correlations as proposed in Ref. [25[ . For the 
vertex function we take the model 



S(q) = 



1 



q 



1 + 



(72) 



where the charge correlation length ^ = 1/ki deter- 
mines the radius of a "correlation hole" . Taking into 
account that £ c /j ~ a/S [25| . we can use the relation 
^ c h = 1/(28) in numerical computations. This gives 
£ c h ~ 10 for the underdoped case (8 — 0.05) and £ c h ~ 2 
for the overdoped case (8 — 0.25). We assume a strong 
EPI g ep = 5t = 2.0 eV and take uj =0.1t= 0.04 eV. 

In computations we use the following parameters for 
the model Q: U = A pd = 8t, t' = -0.2 t" = 
0.10 1. As an energy unit we use t = 0.4 eV. The ex- 
change interaction is described by the function J(q) = 



8 




FIG. 1: (Color online) Electron dispersion in the MFA 
£2(k) along the symmetry directions F(0, 0) — ¥ M(n,n) — Y 
X(ir,Q) -> r(0,0) and X(tt,0) -> Y{0,tv) for <5 = 0.05 (red 
solid line), 0.10 (blue dashed line), and 0.25 (black dash- 
dotted line). Fermi energy for hole doping is at uj = 0. 



2 J (cosq x + cosq y ) with J = 0.4i. For the CI energy for 
the n.n. holes we take V\ = 0.44i and u c = 2.5t. The 
electronic spectrum in the normal state is calculated at 
T = 0.02t ~ 100 K. In computations the grid of 64 x 64 
(k x ,k y ) points and up to 1200 imaginary frequencies oj n 
were used. 



B. Electronic spectrum in the normal state 

At first we consider results in the MFA for the elec- 
tronic spectrum (|25[) . The doping dependence of the 
electron dispersion for the two- hole subband (k) along 
the symmetry directions in the 2D Brillouin zone (BZ) 
is shown in Fig. Q] For small doping, 8 = 0.05, the en- 
ergy at the M(ir, n) and T(0, 0) points are nearly equal 
as in the AF long-range order state. Only small hole-like 
FS pockets close to the (±7r/2,±7r/2) points emerge at 
this doping. With increasing doping, the AF correlation 
length decreases that results in increasing of the electron 
energy at the M(tt, 7r) point and at some critical dop- 
ing 6 ~ 0.12 a large FS appears. At the same time, the 
rcnormalizcd two-hole subband width increases with dop- 
ing from W « 2t at 5 = 0.05 to W « 3t at S = 0.25, 
which, however, remains less than the "bare" Hubbard 
bandwidth W = 4t (1 + 5) where short-range AF correla- 
tions are disregarded. Note that in the dynamical mean 
field theory (DMFT) this narrowing of the subbands due 
to the short-range AF correlations is missed 0, 
while they are taken into account partly in the cluster 
DMFT However, as shown in the DMFT the self- 

energy contribution strongly renormalizes the electronic 
spectrum found in the MFA. 

To consider the self-energy effects in the electronic 
spectrum a strong coupling approximation (SCA) should 
be considered by a self-consistent solution of the system 
of equations for the normal GF (|4U)) and the self-energy 
(fBT)|) . In Ref. [27j a detailed investigation of the normal 



FIG. 2: Spectral function in the SCA along the symmetry 
directions T(0,0) -> M(7T,7r) -> X(ir,0) T(0,0) for hole 
concentration S = 0.10. 




r m x r 

FIG. 3: (Color online) Electron dispersion curves in the 
SCA along the symmetry directions F(0, 0) — > M(n, ir) — > 
X(w, 0) -> r(0, 0) for hole concentration 6 = 0.10. 



state electronic spectrum for the conventional Hubbard 
model in SCA was performed. Therefore, here we present 
only the results of the electronic spectrum computation 
for the model (p} which are important for further studies 
of superconductivity in the model. The spectral func- 
tions (|55p along the symmetry directions are presented in 
Figs. H and H for S = 0.10 and 5 = 0.25, respectively. The 
dispersion curves given by the maximum of the spectral 
function (I55j) at the same doping are displayed in Figs. [3] 
and El 

In comparison with the MFA in Fig. [T] a rather flat 
energy dispersion is found with QP peaks at the FS. In 
general, strong increase of the dispersion and intensity 
of the QP peaks is observed in the overdoped region in 
comparison with the underdoped region. This is in agree- 
ment with our detailed studies of temperature and doping 
dependence of the self-energy (|BT))) and spectral function 
([55")) in [27| which have proved strong influence of AF 
spin-correlations on the spectra. 

To estimate the coupling constant A(q) in the two-hole 
subband, we calculated the renormalization parameter 
Z(q) ([59]) at the Fermi energy, 

Z(q) = Z(q,w = 0) = l + A(q) 
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FIG. 4: The same as in Fig.[2]for hole concentration 8 = 0.25. 




m x r 

FIG. 5: The same as in Fig.[3]for hole concentration 8 = 0.25. 




FIG. 6: (Color online) Doping dependence of the renor- 
malization parameter Z(q) along the symmetry directions 
r(0,0) -5- M(7r,7r) -> X(w,0) ->■ r(0,0) at T w 140 K 
for 8 — 0.05 (red solid line), 8 = 0.10 (blue dashed line), 
8 = 0.15 (pink squares), 8 — 0.25 (black dash-dotted line), 
and 8 = 0.35 (black diamonds). 



= l-[dReS(q,w)/dw]| w=0 



(73) 



The doping dependence of Z(q) is shown in Fig. [SI It 
weakly depends on S in the underdoped case for S < 0.15 
but sharply decreases in the overdoped case for 5 > 0.25. 
The temperature dependence of Z(q) presented in Fig. [7] 
is weak at temperatures lower than the characteristic en- 
ergy of spin fluctuations u> s ~ J. The EPI gives a small 




FIG. 7: (Color online) Temperature dependence of the renor- 
malization parameter Z(q) for 8 = 0.05 at T ~ 140 K 
(red solid line), T w 580 K (black dash-dotted line), and 
T « 1100 K (black squares). Blue dashed line shows Z(q) for 
8 — 0.05 caused only by the spin-fluctuation contribution. 




FIG. 8: (Color online) Doping dependence of the parameter 
X(q). Notation are the same as in Fig. [6] 



contribution to the coupling constant as follows from the 
comparison of Z(q) induced by both spin-fluctuations 
and EPI contributions (red solid line in Fig. [7]) with 
the contribution caused only by spin-fluctuations (blue 
dashed line). The renormalization parameter X(q) (|rj0"|) 
at the Fermi energy 



X(q) = X(q,w = 0) = Re£(q,0), 



(74) 



which determines the shift of the dispersion curve is plot- 
ted along the symmetry directions in Fig. |SJ X (q) de- 
creases with doping in the underdoped region as Z(q), 
but in the overdoped region reveals an irregular behavior 
and becomes small at large doping as Z(q). These re- 
sults demonstrate that at large doping both the electron 
interaction with spin-fluctuations and the EPI become 
weak. 



C. Superconducting gap and T c 

For a comparison of various contributions to the su- 
perconducting gap equation (|56p . as the first step, we 
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consider a weak-coupling approximation (WCA). In the 
WCA, the interaction (l4"4l is approximated by its value 
close to the Fermi energy, \to,z\ ~ 0. Then integration 
over 17 of the dynamical susceptibility in (|4"4"]) yields 



screened CI (l66l we have 



Vc[k) = — > — — . 

N ' q + k 
q 



(83) 



d!7 Imx(q, 11) 

j — z — Q 



d$7 Im%(q, 17) 



= -x(q), (75) 



where x(q) = Rex(q, 17 = 0) is the static susceptibility. 
In the WCA the self-energy contribution in the normal- 
state GF ([49]) is neglected that results in the BCS-type 
equation for the gap function at the Fermi energy y(k) = 
cry>2,a(k,w = 0): 
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w4»-«^tanh^{,( k - q) 



-V(k - q) + [(l/4)|i(q)| 2 + |V(k - q)| 2 ]x c/ (k - q) 
+ | S (k - q)| 2 Xp/l (k - q) 0(w„ - |e 2 (q)|) 
-|i(q)| 2 Xs /(k-q)^ s -| £2 (q)|)}, (76) 

where S(q) = [e|(q) + ^(q)! 2 ] 1 ^ 2 . Whereas for the ex- 
change interaction and CI there are no retardation effects 
and the pairing occurs for all electrons in the two-hole 
subband, the EPI and spin-fluctuation contributions are 
restricted to the range of energies ±uj an d ±w s , respec- 
tively, near the FS, as determined by the (9-functions. 

To estimate various contributions in the gap equation 
(|76j) we consider a model d-wave gap function, ip(k) = 
(A/2) r/(k) where ?7(k) = (cos k x — cos k y ). Then the gap 
equation can be written in the form: 

i = - y [i - 6( q )] 2 ^1! tanh !k j j _ v c 



+V cf + (1/4) |t(q)| 2 x c/ + V ep 9{u Q - |e 2 (q)|) 
-|i(q)| 2 Xs/ 0( Ws -|e 2 (q)|)}. (77) 

In this equation only I — 2 components of the static 
susceptibility and CI give contributions 



V c 



V 



cf 



j^£>(k)cosk, (78) 

k 

l^|l/(k)| 2 Xc/(k) cosfc,, (79) 



Xcf 



v„ 



Xsf 



Sep 

N 



y]5(k) cos fcs, 

k 



(80) 
(81) 
(82) 



Computation yields the following parameters for the n.n. 
intersite CI JB5J: K = Vi = 0.44t w 0.18 eV. For the 



where F c (k) = 0.12 1 (0.28 i) w 0.05 (0.11) eV for 
K = 1 (0.2) , respectively (see Table [TTJ) . Note, that 
the projected CI ([55]) is much smaller than the CI en- 
ergy V c0 (k) = (u c /iV)2]q([l/(? + «)]< In particular, 
K(k)/K (k) = 0.15 (0.24) for k = 1, (0.2), respectively. 
In the conventional BCS theory the CI is suppressed by 
retardation effects described by large Bogoliubov-Morel 
logarithm, hx(jx/u!ph)- In the Hubbard model there are no 
retardation effects for the AF exchange interaction but 
a reduction of the CI contribution is due to the d-wave 
pairing. 

To estimate contributions from the charge fluctuations 
we use the static charge susceptibility (l67|) . Applying 
this approximation to the screened CI (66) we get the 
following expression for charge contribution (79): 



1 



■Xcf s ) cosk x, (84) 



where V cf (n) = 0.05 (0.25) t « 0.02 (0.1) eV for k = 
1 (0.2), respectively. This contribution is smaller than the 
CI repulsion V c (j53"f and in our approximation the d-wave 
pairing induced by the screened CI in the second order 
V c f is destroyed by CI repulsion in the hrst order V c as 
was pointed out in Ref. |19| • The charge fluctuation con- 
tribution from the n.n. intersite CI (j6"5"j) is even smaller, 
V£p ~ 4 • 10~ 3 1 . The contribution from the charge fluc- 
tuations (|80l) calculated for the static susceptibility (|67|) 
is also small: Xc f {8) = (1/t) 0.15-10- 2 (1.3-10~ 2 ) for the 
hole concentrations S = 0.05 (0.10), respectively. For the 
averaged over the BZ vertex |t(q)| 2 = (1/N) ^ q |£(q)| 2 ^ 

4i 2 this contribution is equal to |£(q)| 2 Xcf ^ 0.02 eV and 
can be neglected. 



The EPI contribution (|8 1 1) is given by 



t> _ Sep \ -< £,ch 



cos k x = g ep S d (£ch), (85) 



where S d (£ ch ) = 0.154(0.393) for £, ch = 2(10), respec- 
tively. Thus, even for a strong EPI coupling g ep — 
5i = 2 eV we obtain a moderate contribution from the 
EPI for the d-wave pairing: V ep {£, ch ) = 0.76 1 (1.96 1) w 
0.3 (0.8) eV for £ ch = 2 (10), respectively. The EPI con- 
tribution to the s-wave pairing is given by the I = 
component S = (1/N) £ q S(q) = 0.31 (0.57) for £ ch = 
2 (10), respectively. The ratio of the d-wave S d and the 
s-wave S s — So components of the EPI matrix elements 
is equal to (S d /S ) = 0.43 (0.60) for £ o/l = 2 (10), re- 
spectively. This shows that at small hole concentrations 
S (large charge correlation lengths £ c /j = 1/2(5) the EPI 
for the both components are comparable, while for the 
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TABLE II: CI parameters V c , V c o, V c f, and EPI parameter 
V ep for several values of the CI screening constant k — 45 and 
the charge correlation length £ c f — 1/25 for EPI related to 
hole concentration 5. 



5 


K 




V c /t 


VcO It 


Vcf /t 


Vep It 


0.05 


0.2 


10 


0.28 


1.18 


0.25 


1.96 


0.10 


0.4 


5 


0.22 


1.05 


0.26 


1.4 


0.25 


1 


2 


0.12 


0.80 


0.05 


0.76 



overdoped case the d-wave component Ss, becomes con- 
siderably smaller than the s-wave component in agree- 
ment with the results of Ref . (25[ . 

The spin-fluctuation contribution Xsf calculated for 
the model Xs/(q) in Eq. ([68]) for several values of the 
AF correlation length £ is given in Table HI Using the 
averaged over BZ vertex |*(q)| 2 — 4t 2 we can estimate 
an effective spin-fluctuation coupling constant as g s f ~ 
-At 2 Xs } = 5.3(2.4)* w 2(1) eV for 5 = 0.05(0.25), 
respectively. Thus, the spin-fluctuation contribution to 
the pairing in Eq. (1771) appears to be the largest. Note, 
that g s f is close to the spin-fluctuation coupling constant 
U « 1.6 eV found in Ref. [ll| from ARPES data. 

In Table [TT] we present the coupling parameters in the 
equation for T c (ITTj) . In the MFA the pairing can be in- 
duced by the AF exchange interaction J = OAt = 0.16 eV 
which is comparable with the repulsion caused by the 
screened CI: V c = (0.05 — 0.11) eV or even smaller than 
the n.n. hole CI V\ — 0.175 eV. Therefore, the supercon- 
ducting pairing in the MFA for the t-J model (in partic- 
ular, the RVB state is strongly suppressed (or even 
destroyed) by the intersite Coulomb repulsion. 

To calculate doping dependence of T c we solve Eq. ([77)) 
by taking into account the exchange interaction J, the 
Coulomb repulsion V c , and the contributions from the 
self-energy in the WCA: V ep , Xsf, and V c f, neglect- 
ing the small contribution Xcf- Results of the calcu- 
lation is shown in Fig. [HI The highest T c m 0.2 1 is 
found when all the contributions are taken into account. 
The spin-fluctuation pairing results in superconducting 
T c s/ w 0.04 1 much larger than T c ep w O.Olt mediated by 
the EPI. For the k-independent EPI (5(k) = 1) as in the 
Holstein model, the d-wave pairing is absent. The doping 
dependence of T c is qualitatively agree with experiments 
in cuprates but its value is an order of magnitude higher. 

The high values for T c found in the WCA are explained 
by neglecting the reduction of the quasiparticle weight 
caused by the self-energy effects in the gap equation (1551) . 
It is convenient to write the gap equation in the form: 



T, 



¥>(k,0 = ^££ (j(k-q)-V c (k-q) 

q m 

-Vs/(q, k - q, uj n - u m ) + V ep (k - q,w„ - u) m )} 
[1 -6(q)] 2 <p(q,w m ) 




FIG. 9: (Color online) T c (5) in the WCA induced by all in- 
teractions (red solid line) and only by the spin-fluctuation 
contribution Xs/( Drue dashed line) or only by the EPI V ep 
(black dash-dotted line). 



For V c (k — q) we take the screened CI ([6"6"[) . Since 
the charge-fluctuations gives a much weaker contribu- 
tion than the spin-fluctuation and electron-phonon in- 
teractions (see Table [I] and Table HI)) . we neglect the term 
[|V(k - q)| 2 + |t(q)| 2 /4] x c /(k - q, v n ) in the interaction 
function([64| . Contributions induced by spin- fluctuations 
and the EPI are described by the functions 

F s/ (q,k-q,^) = |i(q)| 2 Xs /(k-q)^/K), (87) 



V ep (k - q,w„) = g ep 



5(k - q) 



where the spectral function for spin fluctuations reads: 

tanh(a; oj s /2T) 



F sf {uj v ) 



1 



2xdx 
x 2 + (w„/w s ) 5 



1 



(89) 



[u m Z(q, Lu m )] 2 + [e 2 (q) + X(q, u m )] : 



(86) 



To calculate T c and to find out the energy- and k- 
dependence of the gap <p(k,u), Eq. ([86]) was solved by 
a direct diagonalization in (k x ,k y , cj„)-space. Since the 
largest contribution in Eq. (|86|) comes from energies close 
to the FS, we have used the renormalization parameters 
at the Fermi energy Z(q) ([73]) and X(q) ([74]) instead 
of the energy dependent ones. The results for T c (5) is 
shown in Fig. [TO] The highest T c - 0.021* - 100 K is 
found when all the contributions are taken into account, 
though pairing induced only by spin-fluctuations also re- 
sults in high T c s/ - 0.014* - 65 K. The d-wave pairing 
induced only by the EPI is rather weak and does not dis- 
played in Fig. 1101 The value of T c is reduced by an order 
of magnitude in comparison with the WCA in Fig. [9] due 
to a suppression of the QP weight by the factor [1/Z(q)] 2 . 
The maximum value of T c is found at lower value of dop- 
ing 5 pt ~ 0.12 than in experiments, S^f = 0.16. 

The k-dependence of the gap function <p(k, uj ~ 0) 
at doping 5 = 0.13 for (0 < k x , k y < 2tt) is plotted 
in Fig. 1111 The gap reveals a distinct d-wave symmetry 
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0.03 




FIG. 10: (Color online) T c (8) induced by all interactions in 
Eq. ()86[) (red solid line) and only by spin-fluctuation contri- 
bution V s f (blue dashed line). 




2 k 4 6 

FIG. 11: (Color online) 2D plot of the SC gap <p(\t,u) ~ 0). 



with maximum values in the vicinity of the FS. As shown 
in Fig. [T2J its angle dependence on the FS is close to 
the model d-wave dependence tpd(6) — cos 20. Energy 
dependence (in units of t ) of the gap function ip(k, oj), 
the real and imaginary parts, is presented in Fig. [T3] at 
k fa (0,7r/2) and 5 = 0.13. Since the gap function was 
obtained as a solution of the linear equation at T = T c the 
value of the gap is given in arbitrary units. The energy 
variation of the gap occurs in the region of uj < 0.4t, 
of the order of the characteristic spin-fluctuation energy 
uj s = J = OAt. 

Generally, the results obtained in the SCA are in a 
qualitative agreement with experiments in the cuprate 
superconductors. They also demonstrate an important 
role of the self-energy effects in the normal and super- 
conducting states in comparison with the MFA. 




40 20 

angle 



FIG. 12: (Color online) Angle dependence of the SC gap <p(8) 
on the FS (blue bold line) in comparison with the model d- 
wave dependence ifid{9) = cos 28 (red dashed lines). 
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FIG. 13: Energy dependence of the real, Re </j(k,w), and 
imaginary, Im ip(k.,uj), parts of the SC gap in arbitrary units. 



D. Comparison with previous theoretical studies 

As briefly discussed in Sec. HI various methods have 
been used in theoretical studies of superconductivity in 
the Hubbard model. Here we would like to emphasize 
our results in comparison with previous investigations of 
the problem. 

At first we refer to results obtained in the weak or 
intermediate correlation limit. In particular, using the 
two-particle self-consistent non-perturbative approach 
(see, e.g. Refs. and the fluctuation-exchange 

(FLEX) approximation (see, e.g., Refs. [5l|, HH and re- 
views [7], |53j), a system of equations was derived within 
the Fermi-liquid model to study self-consistently single- 
electron GFs and the spin and charge dynamical suscep- 
tibility. Within the FLEX approximation, the supercon- 
ducting d-wetve pairing was found in a narrow range of 
doping, very close to the AF instability. In our theory 
superconductivity is mediated by a broad spectrum of 
AF spin excitations (paramagnons) which results in the 
doping dependence of T c close to experimentally observed 
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(see Figs. M HDJ)- 

A general problem of the weak CI in the Hubbard 
model has been extensively studied within the RG ap- 
proach (for reviews see Refs. 0, Ell)- The RG stud- 
ies revealed a competition between various type of 
phases driven by electronic instability, such as the spin- 
density wave (SDW), charge-density wave (CDW), ne- 

matic (Pomeranchuk) phase, stripes, superconducting 

pairing, etc. (see reviews [H, |57[ and Refs. (58U62j ]L 
For a weak Hubbard interaction U ~ t in a certain range 
of hopping parameters and doping the d-wave supercon- 
ducting pairing can overcome other instabilities. An im- 
portant role of the intersite Coulomb repulsion in the 
Hubbard model was found in Ref. [22| . as mentioned in 
Sec.|H In our theory we disregarded other instabilities and 
did not study a general phase diagram since it would de- 
mand investigation of a much more complicated system 
of equations which is beyond the scope of the present 
paper. 

To consider cuprate superconductors, the strong corre- 
lation limit should be investigated. In many publications 
the spectrum of electronic excitations in the normal state 
in the Hubbard model was extensively studied. Here we 
refer to numerical simulations for finite clusters (see re- 
views [63T - [65t ). the DMFT (see reviews EE S3)) tne dy- 
namical cluster approximation (DCA) [66l . |67| and the 
cluster DMFT (see, e.g., Refs. |48|,[68j|). More accurate 
results have been obtained within the DCA and cluster 
DMFT methods where short-range AF correlations are 
partially taken into account. As we have pointed out 
in Sec. MI Al and Sec. IIVB1 in our method short-range 
AF correlations are properly taken into consideration in 
the MFA resulting in a large reduction of the effective 
bandwidth W . Consequently, a two-subband state in the 
Hubbard model is found even in the intermediate corre- 
lation limit, U ~ At. The spectral density computed in 
SCA, Figs. [2] and HI are in accord with numerical studies 
for the Hubbard model (see, e.g., Refs. (Hoi. |6(I |69|). 

The most controversial problem is whether the super- 
conductivity can emerge from the repulsion, as discussed 
in Sec. HI Extensive numerical studies for finite clusters 
have revealed a tendency to the d-wave pairing in the 
Hubbard model, though a delicate balance between su- 
perconductivity and other instabilities (AF, SDW, CDW, 
etc.) was found (see, e.g., Refs. [63468] ) . In Ref. [70j). 
using the DCA with the quantum Monte Carlo method, 
the superconducting d-wave pairing and the isotope ef- 
fect similar to observed in cuprates were found for the 
Hubbard-Holstein model. However, in several publica- 
tions an appearance of the long-range superconducting 
order has not been confirmed (see, e.g., Ref. [z3). There- 
fore, analytical studies are desirable to elucidate this 
problem. 

An accurate analytical method is based on the HO 
technique where the HO algebra is implemented rigor- 
ously. The superconducting pairing induced by the kine- 
matic interaction for the HOs was first proposed by Za- 
itsev and Ivanov [72| who studied the two-particle ver- 



tex equation by applying the diagram technique for HOs. 
The momentum-independent s-wave pairing was found 
in the lowest order diagram approximation equivalent 
to the MFA. However, this solution violates the HO 
kinematics and the t-J model should be used to ob- 
tain the d-wave pairing mediated by the AF superex- 
change interaction (see, e.g., Refs. [73l - [76l |L In this re- 
spect we should point out that in many publications 
superconductivity in the t-J model was studied in the 
MFA (see, e. g., R efs. [z| M, Ezl I?!)- As we have 
shown in Sec. IIV C[ the intersite Coulomb repulsion sup- 
presses or even destroy superconductivity induced by the 
AF superexchange interaction in the MFA. In particu- 
lar, in cuprates, a sufficiently strong n.n. hole repul- 
sion Vi = 0.1 - 0.2 [11 may be detrimental for the RVB 
state [H . The same remark refers to the studies of super- 
conductivity in the conventional Hubbard model in the 
MFA (see, e.g., Refs. [H, [Z^-HU). Therefore, consider- 
ation of the spin-fluctuation pairing beyond the MFA is 
essential in description of superconductivity in cuprates 
as discussed in detail in Sec. IIV CI 

In comparison with studies of the intersite Coulomb 
repulsion in the weak correlation limit in Refs. 0, 
in the strong correlation limit the intersite Coulomb re- 
pulsion Vij to some extent is compensated by the nonre- 
tarded superexchange interaction Jy (seeEqs. (|3"Tj) . (|32|) ) 
which is absent in the weak correlation limit. At the 
same time, even for a sufficiently large component V c q 
of the CI Vij, the contribution to the gap equation is 
given by a much smaller d-wave partial harmonic (|83[) 
and therefore is not so detrimental to superconductivity 
in comparison with the conventional s-wave momentum- 
independent pairing. 

Studies of the spin-fluctuation d-wave pairing in the 
presence of the EPI have shown that depending on the 
symmetry, the EPI could enhance or suppress supercon- 



45] 



ducting pairing (see, e.g., Ref. [82j, |83j). In Ref. 
the d-wave pairing induced by both the spin-fluctuations 
and EPI in the model (fTTj) within the FLEX approxima- 
tion was considered. It was revealed that a momentum- 
independent EPI strongly suppresses T c , while the EPI 
with strong forward scattering can enhance T c . In our 
theory the strong spin-fluctuation pairing is induced by 
the kinematic interaction which is absent in the weak cor- 
relation limit as in FLEX approximation, and therefore, 
the EPI plays only a secondary role in the d-wave pair- 
ing. A strong EPI in polaronic effects observed in the 
oxygen isotope effect on the in-plane penetration depth 
in cuprates [841 ] may be irrelevant for the pairing medi- 
ated by the d-wave partial harmonic of the EPI [85j as 
confirmed by a weak isotope effect on T c in the optimally 
doped cuprates. 



V. CONCLUSION 

In the paper the theory of superconducting pairing 
within the extended Hubbard model |T]) in the limit of 
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strong electron correlations is presented. Using the Mori- 
type projection technique we obtained a self-consistent 
system of equations for normal and anomalous (pair) 
GFs and for the self-energy calculated in the NCA. The 
theory is similar to the Migdal-Eliasberg strong-coupling 
approximation. 

We can draw the following conclusions about the mech- 
anism of pairing in the extended Hubbard model. Solu- 
tion of the gap equation in the weak coupling approxima- 
tion ([75)1 shows that for the d-wave pairing the intersite 
Coulomb repulsion gives a small contribution determined 
by I = 2 harmonic of the interaction function V(k — q). 
However, it can be larger than the AF superexchange in- 
teraction J(k — q) , and the RVB-type superconducting 
pairing can be destroyed. 

Pairing induced by charge fluctuations Xc/i(k — q) ap- 
pears to be quite weak (outside the charge-instability 
region). We have found that the d-wave component of 
the EPI, even for the model of strong forward scatter- 
ing [25| and a large fully symmetric s-wave component, 
turned out to be small. The largest contribution to the 
d-wave pairing comes from the electron interaction with 
spin fluctuations induced by strong kinematic interaction 
|i(q)| 2 , so that the EPI plays a secondary role in achiev- 
ing high-T c . 

It is important to point out that the superconduct- 
ing pairing induced by the AF superexchange interaction 
and spin-fluctuation scattering are caused by the kine- 
matic interaction characteristic of systems with strong 
correlations. These mechanisms of superconducting pair- 
ing are absent in the fermionic models (for a discussion, 
see Ref. (86|) and are generic for cuprates. The intersite 
Coulomb repulsion is not strong enough to destroy the 
d-wave pairing mediated by spin fluctuations. Therefore, 
we believe that the magnetic mechanism of superconduct- 
ing pairing in the Hubbard model in the limit of strong 
correlations is a relevant mechanism of high-temperature 
superconductivity in the copper-oxide materials. 



Note added in proof. 
we became aware of references 



When this work was submitted, 
87|-[89[ which consider 



the extended Hubbard model with the intersite Coulomb 
repulsion V. The results of the references [13], show 
that the on-site repulsion U effectively enhances the d- 
wave pairing which survives for large values of V up to 
V ~ U/2 3> J (Ref. 89]). This observation supports our 
model of spin-fluctuation pairing due to the kinematic 
interaction which emerges only in the strong correlation 
limit. As long as the Coulomb repulsion V does not ex- 
ceed the kinematic interaction of the order of the kinetic 
energy, V < 4t, the d-wave pairing may survive. The 
small value of V = J found in Ref. [H| is explained by 
application of the slave-boson representation in the MFA 
which ignores the kinematic interaction. We would like 
to thank A.-M. S. Tremblay for valuable discussion who 
drew our attention to those papers. 
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Appendix A: Pair correlation function in MFA 



Here we calculate the pair correlation function 
(X® 2 Nj) considering an equation of motion for the com- 
mutator GF Lij(t-f) = ((X° 2 (t) I Ni(t')))\,^. The 
equation for the GF can be written as 



(oj - e 2 )L i:j (uj) 



£ a>tZ{((xr'x° r :'\N J )) u 



((Xf' 2 A*' 2 |A,)U, 



(Al) 



where we neglected excitation energy proportional to the 
intraband hopping in comparison with the interband con- 
tribution, -C I e 2 1 — U . The pair correlation func- 
tion is determined by the equation: 



1 r +ca 

{XfN 3 )=- 1 - - 

7T J -00 1 



da) 



exp(— ui/T) 



ImLy(w). (A2) 



The GF Lij{ui) has two poles, one at w = £2 and an- 
other at the energy of a pair excitation given by the 
GFs at the right-hand side in Eq. (|Alj) of the singly 
or the doubly occupied subbands. Let us consider the 
hole doped case, n = 1 + S > 1 when the chemical 
potential crosses the two-hole subband, fi ~ U and 
£1 = — /i ~ —U, e 2 = 2ei + U ~ —U. In this case we 
can neglect the exponentially small contribution of the 
order of exp(—U/T) <C 1 coming from the pole to = e 2 . 
The contribution from the GF of the one-hole subband 
in {£!]>, 

- -Ira((X^X^\N 3 )) u ~ 5 mj {XfX^ }8(u - 2 £l ), 

(A3) 

also gives an exponentially small contribution of the 
order of exp(— 2U/T) -C 1 . Therefore, we can 
take into account only the contribution from the GF 



((X? 2 X% l 2 \Nj)) LU where the pair excitation energy 



\t 2 ^\. Using the approximation — e 2 ) 



1/U 



which neglects retardation effects an integration over uj 
in Eq. (|A3j) gives the following result, 



-{At\yV)a{XfX-f). 



(A4) 



The last formula is obtained in the two-site approxima- 
tion usually applied for the t-J model: m = j, which 
gives Xf 2 Nj = 2Xf 2 . 
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Appendix B: Self-energy 

The normal and anomalous (pair) components of the 
self-energy operator (I35|) are given by the matrices: 



]T tutfl, ImG 22 CT , (wi)Im«B faff - \Bl a ,)U , (B7) 



ckoidLU2 1 — n(u>i) + N(lu2) 



M ija {uj) = (( 



X Jj H] 



([H,Xf}[H,X9°] 



)<g*\ (Bl) 
>)2*\(B2) 



, t: u> — oji — UJ2 

(B8) 



IV a' 



where [X"^,H] and [H, X"' 13 ] are the irreducible parts 
of the commutators determined by Eq. (|17|) . Using equa- 
tions of motion for the HOs as given, e.g., by Eq. (0 we 
obtain multiparticle GFs which are determined by prod- 
ucts of bosonic and fermionic operators. Let us consider, 
in particular, contributions to the two-hole subband self- 
energy given by the kinematic interaction in Eq. ([7]). The 
normal component in Eq. (|B1I) reads, 



Taking into account the definition of the bosonic operator 
© the bosonic GFs in these equations can be written as 

PWl4"'»<" = (l/4)«iV i |iV i )} w ,W 
+«5?|5|)} a , <W + {{Xf\Xf)) u 5 a , s , (B9) 
(B irTrT ,\B m ,)) u = (l/4)«iV i |iV i )} w <W 
-((Sf \&j)) u tW + ((Xr\Xf )) u cW (BIO) 



M. 



22(fc) 



(7 2 I \ r 2(T // Fit 



jaa" I 



IV a' a 1 ' 



1 f°° , e^ 2 
2 , tatw — I dz 



IV cr' a" 
r2a" n t 



1 



27T ./ „ (w - z) 

cr'2/ 



dte 1 



(B3) 



x {Xf B) aa „\B laa ,(i)Xr{t))- 
For the anomalous component in Eq. (|B2[) we have 



\Biaa' X l 



cr'2 I V S"2 



Xu B~ 



IV a'a" 



IV cr cr 



After summation over cr' in (|B7I) for the bosonic GF 
(|B9I) and the normal GF in the paramagnetic state, 
Gly a {uj) = Gff, s {uj), the spin- fluctuation contribution 
can be written in the form: \{S*\S*)) U + {{Xf a \Xf)) u = 
((Si|Sj)) w . Similar summation over cr' in (|F38|) for the 
bosonic GF (|B10j) and the anomalous GF Ffi? (u)) = 
-Fff s (w), results in the equation: -{{S?\Sfi) u Fff a (u) + 
iX?\Xf)) u F$ s { U ) = -iSi\S s )) u Fff a {u,). 

Introducing the q-representation for the GFs and the 
self-energies as defined by Eq. (TT2"j) for the self-energies 
(|B7I) and (IB8|) we obtain the expressions: 



1 -no 



x (A ; T ^^ s H5„ ff -(f)ir(i)>, 

where the bosonic operator Bi aa i = Bf 2 a , is defined by 
Eq. (JSj> . Using the spectral representation for the ther- 
modynamic GFs [32| we introduced in Eqs. (|B3|) . (IB4|) 
the multi-particle time-dependent correlation functions. 
They are calculated in the mode-coupling approximation 
as described by Eqs. (|40|) . (|4T|) . The time-dependent 
single-particle fermionic and bosonic correlation func- 
tions which appear after the two-time decoupling are cal- 
culated self-consistently as e.g., 



x [-(l/7r)ImG 22 (q,z)], 



(Bll) 



$fW(k )W ) = l^ldz^^zl^k-q) 



x [-(l/t)Im^(q lZ )] 



(B12) 



(Xf, a Xf{t)) = I dw'n^e-^' 1 

J — OO 

x [-(l/^ImG^a/) 

/OO 
du>'N(w')e 
-OO 



where the contribution from the kinematic interaction is 
given by the kernel 



iuj t 



\B 



(B5) 



,(B6) 



x [-(l/7T)]Im<<IW|^ W , 

Here Gf^ a {uj') is the GF J3T| for the two- hole subband 
and ((BiCTCT'|-Sj CTCT '))w is the commutator GF for bosonic 
excitations. Integration over time t in Eqs. (|B3|) and (|B4p 
yields 

M 22(fe)^-j _ f f°° duj x duj 2 1 - n(cji) + iV(a;2) 



lr (±) l , . v l<(q)| 2 /• Jo l + jV(0)-n(z) 

W TT J U) — Z — il 

— oo 

x {lm Xs/ (k - q,n) ± (l/4)Im Xc /(k - q,n)}. (B13) 

Here the spin- and charge-susceptibility are defined by 
Eqs. (l4"5j) and (|4^|) . By taking into account contributions 
from the CI and the EPI in Eq.© and using the NCA 
in calculation of the respective time-dependent correla- 
tion functions we obtain the kernel (14"4"|) for the integral 
equations 02} and ([i3"j). 
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